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Abstract 

A  generic  computer  code  based  on  the  iterative  boundary  integral  equation  method  (I-BIEM)  is 
developed  for  simulating  a  variety  of  electrochemical  problems.  In  this  work  we  have  extended  the 
reach  of  the  method  by  developing  a  generalized  program  capable  of  solving  a  wide  range  of  elec¬ 
trodeposition  problems.  The  new  code  accomodates  quite  easily  multi-variable  problems  including 
those  with  curved  boundaries,  and  non-linear  boundary  conditions.  Such  problems  include  anoma¬ 
lous  c.odeposition  of  alloys,  incorporating  the  effects  of  convective  and  diffusive  mass  transport, 
tittle  variant  effects  such  as  would  be  observed  in  extended  growth  calculations  and  pulse  plating: 
tnicrostructural  modeling  with  reference  to  nonisotropic  boundary  conditions  and  crystallographic 
effects.  Some  interesting  results  f  f  real-life  simulations  are  presented. 
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Introduction 


Complex  electrochemical  processes  in  domains  of  nontrivial  geometries  are  generally  difficult  phe¬ 
nomena  to  simulate.  The  advent  of  the  iterative  boundary  integral  equation  method,  1-BIEM, 
developed  by  Cahan  et  al  (1988),  has  made  the  efficient  study  of  such  complicated  processes  a 
possibility  using  a  computer  of  medium  memory  and  storage  capability.  The  flexibility  that  the 
method  allows,  in  terms  of  model  refinement  at  no  significant  increase  in  memory  or  storage  require¬ 
ments,  gives  this  new  technique  a  major  advantage  over  existing  codes  for  studying  electrochemical 
processes. 

Convective-diffusion  processes  are  handled  efficiently  in  a  way  that  minimizes  numerical  dis¬ 
cretization  errors  that  commonly  pollute  the  description  of  the  underlying  physical  process  when 
using  other  numerical  methods. 

Since  all  essential  numerical  operations  are  performed  solely  o*'  prohleo?  boundary,  the 
method  is  an  excellent  tool  tor  modeling  codeposition  and  microstructuring.  This  is  because  the 
computed  variables  are  precisely  those  pertinent  boundary  quantities  required  for  the  solution 
The  iterative  nature  of  the  method  makes  for  an  efficient  handling  of  transient  phenomena.  The 
dynamic  changes  in  geometry  often  occurring  in  electrochemical  processes  as  in  the  formation  of 
dendrites,  nodules  and  other  surfaces  including  fractals  are  neatly  simulated  using  the  I-BIEM. 

Many  of  these  problems  are  governed  by  Poisson-type  equations  rather  than  the  Laplace's 
equation.  In  the  latter  the  right  hand  side  is  zero  whereas  we  use  the  term  Poisson-like  to  describe 
all  other  problems  where  the  right  hand  side  may  be  a  constant,  a  known  function  of  position 
and/or  time  or  a  function  of  the  dependent  variable  it.self. 

In  the  development  of  the  new  generic  code  the  following  steps  have  been  taken: 

1 .  Several  practical  electrochemical  problems  are  described  in  terms  of  fundamental  or  empiri¬ 
cally  determined  relationships; 

2.  Appropriate  differential  and/or  integral  equations  describing  the  system  are  selected 

.3.  Algorithmic  tools  are  then  developed  to  incorporating  these  into  1-BIEM. 

Description  of  Key  Electr  ^  •'mical  Problems 

•  Coupled  Diffusion-Current  i  '  oution/Anamolous  Codeposition 

In  many  cases  simple  solution  ol  the  current  distribution  of  the  problem  without  consideration 
of  other  variables  can  provide  enough  information  to  give  deep  insights  into  the  underlying 
principles  of  a  particular  problem.  (This  is  often  helped  by  the  fact  that,  for  example,  diffusion 
problems  have  solutions  similar  although  not  exactly  those  of  potential  distribution  problems, 
at  least  in  small  dimensions.)  Intermediate  scale  problems  can  often  be  handled  by  the  use 
of  simple  “correction  terms”  based  on  well  recognized  engineering  principles.  In  some  cases, 
such  as  in  alloy  deposition,  this  is  more  often  than  not  the  weak  point  of  a  treatment. 

In  such  ca.ses  a  more  rigorous  solution  involves  the  simultaneous  solution  of  multiple  variables 
which  are  interrelated  at  the  boundaries.  A  rigorous  treatment  of  anomalous  deposition,  for 
example,  will  require  simultaneous  solution  of: 

1.  potential-current  distribution  ; 

2,  concentrations  of  two  independent  species  (Ni  and  Co)  and 


3  local  pH. 


Only  (1)  can  be  treated  using  a  straight  Laplacian  formulation.  The  others  have  to  be  given 
various  levels  of  sophistication  in  Poisson-like  formulations.  For  example 

-  Simple  diffusion  in  a  large  scale  system.  Here  only  the  region  near  the  boundary  is  rele¬ 
vant.  Most  of  the  bulk  liquid  is  relatively  unaltered  by  the  electrodeposition  processes. 

-  Several  of  the  species,  such  as  Ni  and  Co,  can  be  treated  as  independent  species  connected 
only  by  the  boundary  (kinetics  etc.)  equations.  No  cross  terms  in  the  bulk  need  be 
considered.  If  the  pH  is  a  relevant  variable  and/or  the  precipitation  of  films  occurs  at 
the  interface  cross  terms  for  the  bulk  must  now  be  considered. 

•  Composition  Variation 

Since  the  boundary  equations  are  different  for  each  of  the  sub  fields  discussed  above,  the 
problem  of  lateral  variability  of  such  deposits  is  not  necessarily  obvious.  When  compositional 
variation  in  the  thickness  of  the  media  is  important  this  must  be  treated  by  developing  the 
time  dependent  solutions  of  the  sets  of  equations  described  above. 

•  Microstructural  Models 

To  date  I-BIEM  has  been  used  primarily  for  macroscopic  modeling  studies.  Since  the  same 
fundamental  governing  ecpiations  apply  in  the  micro  scale  with  the  exception  that  stochastic 
variation  in  local  parameters  must  be  considered.  In  addition  the  properties  of  the  deposit 
including  crystallographic  orientations,  the  development  of  grain  boundaries  and  the  elfecis 
of  localized  segregations  in  the  bulk  and  at  the  grain  boundaries  must  be  considered. 

•  Induced  Codeposition 

1'he  prttblem  in  above  is  further  complicated  when  minority  species  are  to  be  considered.  Here 
the  local  chemistry  and  or  physics  of  the  situation  must  be  incorporated  into  the  boundary 
equations. 

•  Extended  Crowth 

The  mathematical  formulation  has  been  shown  to  be  stable  even  in  the  presence  of  physically 
unstable  conditions.  The  local  effects  of  real  surface  roughness,  impurity  inclusions  and  other 
similar  profrlems  can  be  treated  as  extensions  of  the  above.  From  a  fundamental  point  of 
view  as  well  as  a  practical  concern  the  interaction  of  the  relevant  variables  is  interesting 
Such  studies  could  for  example  lead  to  a  better  understanding  of  the  formation  of  nodules, 
dendrites  or  other  irregularities  in  printed  circuit  fabrication 

•  (Complex  Geometries 

Up  until  now  primary  emphasis  in  the  development  of  the  l-BIEM  has  been  placed  on  the 
mathematical  and  the  algorithmic  portions  of  the  code.  In  order  to  make  the  system  more 
“user  friendly”  and  thus  applicafrie  to  a  broad  range  of  problems  it  was  desirable  to  devehip 
a  menu  driven  method  of  readily  entering  or  altering  boundary  geometries  and  governing 
equations. 

•  fhirrent  Thieves  and  Auxiliary  Anodes 

1'his  problem  is  readily  handleable  once  the  problem  of  assignment  of  complex  geometries  is 
solved. 
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•  Stresses  in  Films/Convective  Transport 

Stresses  in  films  are  often  a  consequence  of  sub  processes  like  hydrogen  codeposition  in  elec- 
trodeposited  films.  In  such  films  the  bi-harmonic  equation  can  be  used  to  treat  these  effects. 
The  bi-harmonic  equation  should  be  solvable  in  terms  of  the  equivalent  set  of  simultaneous 
Laplace  and  Poisson  equations.  These  same  equations  could  be  used  to  study  the  vorticity 
inherent  in  the  convective  transport  in  plating  baths  caused  by  rotation  of  the  disk. 

Governing  Differential/Integral  Equations 
Preamble 

The  class  of  problems  outlined  above  are  generally  governed  by  linear  second  order  partial  differ¬ 
ential  equations  of  the  Poisson-type  The  material  coefficients  involved  in  some  of  the  cases  are 
functions  of  the  space  coordinate  (heterogeneity).  For  such  problems  the  governing  equation  can 
be  cast,  using  parameter  perturbation  technique  in  a  form  that  is  well  amenaole  treatment  by 
the  I-BIEM.  The  result  will  be  series  of  Poisson  equations  the  solution  of  each  of  which  represeni,.s 
identifiable  components  of  the  material  inhomogeneity. 

There  are  a  number  of  approaches  for  treating  terms  involving  time  derivatives  in  transient 
problems: 

1.  The  most  direct  approach  is  the  use  of  special  time-dependent  fundamental  solutions  in  t  he 
integral  equations  for  the  1-BIEM  code.  Such  solutions  are  different  from  the  usual  logarithmic 
solutions  used  in  integral  equations  derived  for  straight  Poisson-type  equations.  However 
the  use  of  such  complex  fundamental  equations  call  for  equally  tedious  implementation  of 
numerical  integrations  in  the  computer  [jrogram. 

2.  The  unsteady  governing  equations  can  be  transformed  so  that  time  is  in  effect  integrated  out 
A  suitable  transformation  is  the  l.aplace  transform  method.  This  will  convert  the  partial 
differential  equation  into  a  ffemholl z-type  of  equation  The  fundamental  solution  in  this 
case  is  in  the  form  of  a  Bessel  function.  The  main  challenge  in  this  approach  is  tfie  inverse 
transformation  back  to  the  physical  time  domain.  There  exist  in  the  literature  a  number  of 
efficient  numerical  algorithm  (e  g.  Schapery  1962,  Bellman  em  et  al  1966,  Stehfest  1970)  for 
performing  the  inverse  transformation. 

3.  The  terms  involving  time  derivatives  are  replaced  by  their  equivalent  finite  differences.  The 
resulting  equations  are  either  of  the  Poisson  or  Helmholtz  type  and  the  numerical  process 
will  necessarily  involve  performing  area /domain  integrations  in  addition  to  the  calculations 
carried  out  on  the  boundary. 

Mathematical  Formulation 

•  Differential  Equations 

Consider  an  electrochemical  domain  H  consisting  of  the  boundary  P.  In  general  the  process 
can  be  described  by  the  differential  ecpiation: 

V^«|>(X,  t)  ^  ^y-X.t)  (  I  ) 
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where  V*  is  the  Laplacian  operator  in  space  (x),  t  is  time,  4>  can  represent  voltage  (V^), 
concentration  (c)  or  temperature  (T)  while  7  represents  some  forcing  functions  which  takes 
on  the  following  forms  depending  on  the  character  of  the  problem: 


-  Diffusion  Problems  c 


7  =  D 


where  D  is  the  diffusion  coefficient.  In  that  case  eqn(  1)  when  transformed  into  the 
Laplace  domain  takes  on  the  form: 

-  l3^{c  -  ct=o)  =  0 

in  which  0  =  \/'sD,  s  is  the  transform  parameter  and  c  =  cexp(-st)  dt. 

Convective  Problems 


7  =  ~v  •  ^<t> 


where  v  is  the  velocity. 

Heat  Production  (Sources  A:  Sinks) 


where  S  is  the  Dirac  delta  function  while  /,  represents  the  strength  of  the  t-the  source 
or  sink. 


Non-linear  Problems 


r  ^  ri 


The  conditions  associated  with  the  boundary  P  can  in  general  be  written  in  the  form: 

(7rj<3i(x)]  +  -  0  (2) 

where  P  is  a  differential  operator,  <7  is  a  coefficient  and  /  is  some  prescribed  function  of  (t> 
On  a  conductor  0=0  while  on  an  insulator  /  =  0. 

Integral  Equations 


The  first  task  in  a  BIEM  process  is  to  convert  the  governing  partial  differential  eqn(  1)  into 
a  suitable  integral  equation.  This  is  achieved  by  multiplying  the  equation  by  a  function  g. 
integrating  the  ensuing  expression  over  H  and  invoking  the  Green’s  identities.  The  result  of 
such  manipulations  is  the  integral  equation: 

a</>(x)  =  f  {<^(x')^~(x,x')  ^  9(x,x')^*’'’(x')}  (ir(x')  f  f  g{x,x")?'(x")  dn{x")  (3) 
Jr  on  on  Ju 

where  n  is  along  the  unit  outward  normal  to  the  boundary  and  a  is  the  Cauchy  principal 
value  of  the  integration  of  the  (Jreen's  function  singularity  since  the  following  equation  is  to 
be  satisfied  by  g: 

£[(7(x,x')|  =  ^(x,x')  (4) 

where  £  —  for  Poisson-type  problems  and  £  =  -  0^  for  Helmholtz-types 
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I-BIEM 

Coefficient  Matrices 


To  use  the  boundary  integral  equation  method  the  boundary  F  is  subdivided  into  a  finite  number 
of  elements  and  suitable  interpolation  functions  are  chosen  to  represent  the  distribution  of  the 
dependent  variables  on  the  boundary.  For  function  ^  on  F  we  write: 

=  12  (^) 

t 

where  are  the  values  of  at  the  discrete  points  on  the  boundary  and  iV,  are  shape  functions.  A 
similar  expression  can  be  written  for  d4>ldn.  To  perform  boundary  integrations  we  select  a  finite 
number  of  points  to  serve  as  integration  origins  (x  in  eqn  3).  In  most  BIEM  implementations  these 
points  are  located  along  the  boundary  segments  and  usually  fall  on  the  nodes  formed  at  element 
intersections.  In  a  few  implementations  (the  so-called  non-singular  formulation)  the  integration 
points  are  selected  outside  f2.  In  I-BIEM  x  can  be  placed  on  F,  inside  or  outside  of  Q  The 
number  of  independent  integration  points  selected  is  the  same  as  the  number,  Ai,  of  unknown  </> 
and  <t>n  -  d(f>/dn  on  the  boundary. 

When  these  integrations  are  performed  eqn(  3)  becomes: 


where 


M 

r^i 


M 


<hj  -  (x')  dr 

f’,,  -  ^^(x,x')/V,(x')  dF 

c,  =  /  ^(x,x")T{x")  dC(x") 

Jn 


(fi) 


W'hep  the  integration  points  x,  fall  outside  H  then  =  0. 

To  solve  for  the  unknown  4>  and  4),,  at  the  nodes  we  start  by  guessing  any  initial  values,  and 
(t>n  for  these  quantities.  On  a  conductor  (exact  <t>^)  or  on  an  insulator  (exact  )  type  problems 
the  known  quantities  can  be  used  once  and  for  all  in  e<]n(  fi)  and  the  results  absolved  into  the 
coefficients  c, .  For  linear  mixed  boundary  conditions  one  will  normally  eliminate  either  (i>  or  <i>„  for 
the  other.  In  non-linear  mixed  conditions  the  iterative  process  readily  accepts  the  incorporation  of 
any  root-finding  routine  in  the  solution  prore.ss  as  will  become  clear  shortly  When  the  guesses  are 
used  in  the  equation  written  for  the  «-th  node  as  origin  of  integration  (see  eqn  6)  the  result  is  an 
error: 

\f  M 

(7) 

;-i  7-1 

We  then  update  the  guesses  by  assuming  on: 
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•  Flux  Boundaries  (e.g.  Insulators) 

4>i  —  <i>‘i  ^  (S) 

On 

•  Dirichlet  Boundaries  (e  g.  Conductors) 

<f>n.  =  K.  '  A."’ 

•  Mixed  Boundaries  (e  g.  Tafel) 

a..(d>.  -  4>\)  ~h„(4>n.  -  =  Af.  (10) 

where  A  is  in  an  over-relaxation  factor  For  a  segment  with  a  mixed  boundary  condition  eqn(  in) 
is  to  be  solved  with  the  relevant  prescribed  condition  as  shown  by  eqn(  2). 

The  above  process  is  repeated  for  all  points  until  a  suitable  convergence  criterion  is  satisfied. 

Capabilities  of  the  Developed  I-BIEM  Code 

•  Poi,«,oon  Equation.!! 

The  code  handles  the  solution  of  Poisson-type  equations.  This  permits  the  evaluation  of 
the  effect  of  diffusion  convection  and  time  varying  problems.  Problems  other  than  these 
mentioned  above  which  could  be  tackled  with  this  module  would  include; 

Pulse  plating; 

Magnetic  fields  including  detailed  modeling  of  recording  heads; 

Sources  and  sinks  as  might  be  found  in  problems  associated  with  internally  generated 
thermal  fluxes 

•  Tmsplntinn 

In  order  to  deal  with  problems  governed  by  Poisson-type  equations  it  was  necessary  to  develop 
a  generalized  and  optimizeable  scheme  for  subdivision  of  the  domain  (however  complex)  into 
readily  integrable  subdomains.  This  has  been  devised  to  require  a  minimal  intervention  from 
the  user’s  point  of  view. 

•  Generalized  Boundary  Input 

Automatic  definition  of  the  position  and  subdivision  of  the  boundaries  for  particular  problems. 
Provides  a  generalized  CAD-like  data  input  approach  for  a  wide  variety  of  problems. 

Examples 

Conclusions 
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